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1 Introduction 

Complex systems of many different kinds are in the focus of modern research 
[1,2,3,4]. Correlation matrices obtained from data sampling are a key tool 
to study such systems [5,6]. A standard approach in multivariate statistics is 
to use random matrix theory (RMT) to model the correlation matrices [7]. 
Adopting the framework of RMT, we will calculate the marginal probability 
density function for any single eigenvalue to take a given value x (referred 
to as the one-point function Sp(x) for short) for the case of real ((3 = 1) 
and complex (f3 = 2) correlation matrices. In the complex case, an explicit 
expression is known [8] . Real correlation matrices are encountered more fre- 
quently, but unfortunately, closed expressions for their one-point functions 
and related quantities have not been obtained so far. This is because a certain 
integral over the orthogonal group is not available in explicit form. Sophis- 
ticated power series techniques have been developed in order to tackle this 
problem [7]. However, the resulting expressions suffer from the drawback 
that a resummation of the infinite series has not been possible so far. For 
correlation matrices of large dimension, asymptotic results were derived in 
[9]. Recently some new results for the one-point function and the two-point 
correlation function in the asymptotic regime have been found [10]. 

Here we provide exact results for the one-point functions of real and 
complex correlation matrices. We use an alternative approach to circumvent 
the problems mentioned above. Our approach relies on the supersymmetry 
method [11] - nowadays a standard tool for RMT applications in physics [12]. 
We derive an exact expression for the one-point function of real correlation 
matrices as a twofold integral. We also rederive the known result for the one- 
point function of complex correlation matrices. In Rcf. [13], we presented our 
main results to make them available for applications. Here, we give the full 
derivation of our results in a form that addresses not only physicists and 
practitioners, but also the mathematics and statistics community. 

The article is organized as follows. We formulate the problem and intro- 
duce our notation in Sec. 2. In Sec. 3 we apply the supersymmetry technique, 
pursuing in parallel two different approaches put forward in the literature, 
namely the generalized Hubbard-Stratonovich transformation [14, 15] and the 
superbosonization formula [16,17]. These two approaches are equivalent [18], 
but they yield different forms of the final expressions, each having their own 
advantages and disadvantages. We explicitly calculate the one-point function 
in Sec. 4. In Sec. 5 we numerically integrate our formula to compare with 
Monte-Carlo simulations. We summarize and conclude in Sec. 6. 



2 Formulating the Problem 

In Sec. 2.1 we define the ensemble of random matrices to be considered in 
this paper. In Sec. 2.2 we introduce a generating function for the one-point 
function. This generating function will serve as the starting point for the 
supersymmetry method. 
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2.1 Ensemble of Wishart correlation matrices and one-point function 

We briefly sketch the RMT approach to correlation matrices as set up in 
Ref. [7]. We consider real and complex Wishart correlation matrices. The 
building block for these are rectangular p x n matrices which we denote by 
W = [Wjk], with j = l,...,p and k = 1, . . . ,n. The p rows can be viewed 
as the model time series of length n. The p x p matrix WW^ is the model 
correlation matrix, also referred to as the Wishart correlation matrix. We 
always assume p < n. The entries of W are either real or complex random 
variables. These two cases are labeled by the Dyson index (3 taking the value 
= 1 for real entries (Wjk € M) and /3 = 2 for complex entries (Wjk € C). 
For the joint probability distribution of the entries of W one chooses the 
multivariate Gaussian distribution 

Pp(W,C) = D f3 e X p(-^tiWiC- 1 w), (1) 

where C is the empirical, i.e. given correlation matrix. By construction, the 
ensemble-averaged Wishart correlation matrix WW^ /n equals C, the empir- 
ical one. The Gaussian assumption (1) is justified in most if not all situations 
of interest. The full measure is Pp(W, C)d[W] where 

( V n 



d[W] = < 



n n dw i k f ° r p =i > 



'ft (2) 
j[ n dReWj k dlmWj k for (3 = 2, 



is the corresponding volume element. This measure fulfills the invariance 
condition 

Pp(W,C)d[W] = Pp(UW, UCU r )d[UW] (3) 

for an arbitrary orthogonal ((3 — 1) or unitary (ft = 2) p x p matrix U . 
Since the domain of W (M pxn for j3 = 1, C pxn for /3 = 2) is invariant under 
the transformation W t— > UW , we may replace C by the diagonal matrix of 
its eigenvalues as long as invariant quantities such as the one-point function 
(see below) are studied. Thus, we set C = A where A = diag(yli, . . . , A p ) is 
the diagonal matrix containing the eigenvalues. By the definition of C as a 
correlation matrix we have Aj > 0. The constant Dp in Eq. (1) ensures the 
normalization of Pfj(W, C)d[W] to unity and is given by 

D P = ((2n/0) p detAy n0/2 . (4) 

The set of random matrices WW^ with the entries of W distributed according 
to Eq. (1) is referred to as the ensemble of Wishart correlation matrices 
(sometimes also as the correlated Wishart ensemble). We mention in passing 
that for the choice A = t p , where l p denotes the p x p unit matrix, the 
ensemble defined by Eq. (1) is equivalent to the so-called Gaussian chiral 
random matrix ensemble, which serves as a model for the universal eigenvalue 
statistics of the Dirac operator in Quantum Chromodynamics [19]. 
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The one-point function for the eigenvalues Aj of WW^ is defined by 

r i p 

Sp(x) := / d[W}Pp(W,A)-J2Hx-X J ) 

= — lim Im / d\W]Pa(W, A) tr-. . * p , (5) 

7rpe^o+ J (x- ie)l p - WWt v ; 

where in the second line we have passed to an expression involving the resol- 
vent. By the definition (5), the one-point function is a function of x which 
parametrically depends on the empirical eigenvalues A\, . . . , A p . We drop the 
dependence on A in writing Sp(x). 

Having defined the object of interest, we briefly comment on why it is diffi- 
cult to handle the real case /3 = 1 by the traditional techniques of multivariate 
analysis. In the standard approach to calculating the one-point function (5) 
one makes a singular- value decomposition 

W = UwV, (6) 

where U G 0(p), V G 0(n) for ft = 1 and U € U(p), V € U(n) for /3 = 2. 
The px n matrix ?z; contains the singular values Wj € R (J = 1, . . . ,p) of W. 
For the matrix WW' this decomposition yields 

t^M/t = ^2^1 with w 2 = diag(w 2 , . . . , w 2 p ). (7) 

While the substitution W h- > C/M^ leaves the resolvent in Eq. (5) invariant, 
the diagonalizing matrix {/ does not drop out of the probability distribution 
function P/3(W, A) containing the matrix A ^ l p . The decomposition thus 
leads to the group integral 

$fi(A,w 2 )= fexp(- ^trrf A- 1 Uw 2 ^dn{U). (8) 

For j3 = 2 this is the celebrated Harish-Chandra-Itzykson-Zuber integral [20, 
21] for which there exists an explicit expression. On the other hand, for /3 = 1 
no simple expression is known. The only expression available [7] is an infinite 
series in terms of zonal polynomials or equivalently, Jack polynomials, which 
in turn are only known recursively. A resummation of the infinite series has 
not been possible so far. In the present work we circumvent this problem by 
using a supersymmetry approach. This allows us to derive an expression for 
the one-point function as a twofold integral. 



2.2 Generating function 

The starting point for the supersymmetry approach is the generating function 
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where Xq, X\ are complex variables, Xq ^ M + . The one-point function can be 
computed from it by taking a derivative: 



ax i 



lim (Zp(x — ie, x\) — Zp(x + ie, xi)\ . (10) 



Note that the generating function Z^(xq 1 x\) equals unity at xq = x\. In the 
following we derive simple and computationally useful expressions for it by 
applying the supersymmetry technique. 



3 Passing to superspace 

There are several ways to express the generating function as an integral 
over a suitable superspace. Of particular prominence are the generalized 
Hubbard-Stratonovich transformation put forward in Ref. [14, 15] and the 
superbosonization formula derived in Ref. [17]. Superbosonization [22] was 
first proposed in a field theoretical context. It was then explored how the su- 
persymmetry method can be extended to arbitrary invariant random matrix 
ensembles. In the unitary case, this problem was solved [14] by introducing 
the generalized Hubbard-Stratonovich transformation. In Ref. [17], rigorous 
superbosonization was developed for all classical Lie symmetries (unitary, 
orthogonal, symplectic). The approach of Ref. [14] was then completed by 
transcribing it to the orthogonal and symplectic cases [15]. The equivalence of 
the superbosonization of Ref. [17] and the generalized Hubbard-Stratonovich 
transformation of Refs. [14,15] was demonstrated in [18]. 

In Sec. 3.1, we write the ratio of determinants in the generating func- 
tion as a Gaussian integral over a rectangular supermatrix. Then we carry 
out the ensemble average. To the reader not experienced with anticommut- 
ing variables, we recommend the introductory parts of Refs. [11,23,24] and 
the book by Berezin [25]. In Sec. 3.2, we use a duality between dyadic ordi- 
nary matrices and dyadic supermatrices to express the result of the ensem- 
ble average as a supermatrix integral. After analyzing certain symmetries of 
this dyadic supermatrix we replace it by a supermatrix of the same symme- 
tries but independent matrix elements by means of the generalized Hubbard- 
Stratonovich transformation (Sec. 3.3), and alternatively with the help of the 
superbosonization formula (Sec. 3.4). 



3.1 Ensemble average 

The determinant in the denominator of Eq. (9) can be expressed as a Gaus- 
sian integral over a vector comprising ordinary commuting variables. The 
determinant in the numerator can be expressed as a Gaussian integral over a 
vector with anticommuting entries [23, 11] . By combining both expressions we 
obtain a representation for the ratio of determinants in Eq. (9) as a Gaussian 
integral over a rectangular supermatrix A comprising both vectors: 

det ( Xl l p - WW^) f d[A] ^ /i| {xA ^ A _ AiwwiA) \ (n) 



det(x l p - WW*) J L J \ 2 J 
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where the diagonal matrix X is given by X = diag(a;o, Xq, X\, X\) for f3 = 1 
and X = diag(xo, x{) for j3 = 2. We take Imxo > in order for the Gaussian 
integral to converge. The rectangular supermatrix 



A= [u a ,V a , C,Ca] 1<a< _, ^ 



A 



[ z a (a] i< a <„ > 



At 



Ub 
Vb 

Cb 

-a 

Zb 

S b i i< b <p 



for 0=1, 



l<b<p 



(12) 



for /? 



is p x (2//3|2//3) dimensional. Here € K and Zj € C are ordinary real 

or complex variables while Cj>C| are anticommuting variables, also referred 
to as Grassmann variables. We denote by z* the complex conjugate of Zj. In 
Eq. (11), d[A] denotes the following product: 



d[A] 



(2tt)-p fl du j dv 3 d q d < 3 for P = 1 , 



j'=i 



7r p FT e£Re 2j dim Zj d^+d^. for /3 = 2. 

3 = 1 



(13) 



It should be mentioned that, in this context, one often writes dC = for 
a Grassmann variable £. However, the transformation law d(tQ = t~ 1 d( ! for 
t e C shows that d( = d/d( really is a partial derivative, not a differential! 

By inserting the representation (11) into the generating function (9) and 
changing the order of doing the integrals we find 

Zf}(x ,xi) = I d[A}exp ( ^stvXA^A) 



xD, 3 J d[W] exp ^-^trW^A" 1 + iAA*)W 
j d[A] exp strXA^Aj det- ,l,3/2 (l p + iAAU) . (14) 



In the last step, we performed the Gaussian integral over W. 



3.2 Duality between ordinary and superspace 

We now rewrite the determinant in Eq. (14) as a superdeterminant. This is 
possible due to duality relations between ordinary spaces and superspaces, see 
Refs. [26,14,15]. In the present context the duality amounts to the relation 
[14,15] 

det (l p + iAA^A) = sdet (l 4//3 + iA^AA) . (15) 

We notice that the determinant is a polynomial while the superdeterminant 
is in principle a rational function. The relation originates from the identity 
tr (A A* A) = str (A* A A) and a Taylor expansion in the Grassmann variables, 
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which is always a finite sum. The supermatrix A^AA has dimension 4x4 
and 2 x 2 for f3 — 1 and f3 = 2, respectively. On the other hand, the original 
matrix AA' is p x p dimensional. This dimensional reduction is the crucial 
advantage of the supersymmetry method. 

For present use, we take a look at the symmetry properties of the super- 
matrix A^AA. We see that in both cases (/3 = 1, 2) the left upper block (a.k.a. 
the boson-boson block) of A^AA is a Hermitian matrix. This observation will 
constrain some of the matrix blocks appearing below. What about the com- 
plex linear symmetries (i.e. those not involving complex conjugation)? For 
ft = 2 there are no such symmetries, but for j3 = 1 we have 



(16) 



reflecting the fact that the related pxp matrix A 1 / 2 AA^ A 1 / 2 is symmetric. 

Our aim now is to replace the supermatrix A^AA by a supermatrix a 
with independent matrix elements. We have two approaches at our disposal: 
the generalized Hubbard-Stratonovich transformation [14,15], and the super- 
bosonization formula derived in Refs. [16,17]. 



{A^AAf = S T A^AAS with S = 



10 
10 
-1 
10 



3.3 Generalized Hubbard-Stratonovich transformation 

We proceed by introducing a (super-)Fourier representation of the required 
power of the superdeterminant function on the right-hand side of Eq. (15): 



(17) 



sdet-"' 3/2 (l 4/ , 3 +\A*AA) = J d[g}I p (g)exp (-^stiA^AAi 

The Fourier transform Ip(g) is 

Ip(e) = J d[a]sdet~ nfi/2 (l i//3 +kj)exp strcr^ . (18) 

In order for this integral representation to be formally consistent, the two 
supermatrices a and g have to share the complex linear symmetries of A^ AA. 
Hence a and g are supermatrices of dimension 4x4 for (3 — 1 and 2x2 for 
/3 = 2, and for /3 — 1 the complex linear constraint (16) is imposed. 
We write the supermatrix a as 



Co X 



(19) 



where the entries are 2//3 dimensional square matrices. The matrix cto (akin 
to the boson-boson block of A'AA) is Hermitian, and o\ is a real scalar. 
(The reason for putting the imaginary unit in front of o\ will become clear 
presently.) The off-diagonal blocks \ an d X contain all anticommuting vari- 
ables of a. For /? = 2 the diagonal block <7q is simply a real number whereas 
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X and x = X* are two Grassmann variables. For (3 — 1 the diagonal blocks 
Co and <J\\i are real symmetric 2x2 matrices. The off-diagonal blocks for 
(3=1 have the following structure: 



X 



X 



(20) 



where r; and £ denote Grassmann variables. With this choice a satisfies the 
constraint (16). The supermatrix g, similar to a, is chosen as 



Q = 



g uj 
Co IQlt-i/p 



(21) 



In a self-evident way, g is divided into blocks having the same symmetries as 
the blocks of a in Eq. (19). 

The super-integration measure for a in Eq. (17) is flat and reads 



d[a\ 



(2tt) 2 d<r 0aa d(T ab da obb da x d v d v *d^d^ 
(27r) _1 d(To dai d x d x * 



for (3 = 1, 
for (3 = 2, 



(22) 



where ao aa and cofcb are the diagonal elements and <7o a 6 is the off-diagonal 
element of the real symmetric matrix <7o for (3 = 1. The measure d[g] is 
defined by an identical expression. 

We now insert the representation (17) into Eq. (14). Our ensemble-averaged 
generating function then becomes a supermatrix integral: 



d[g] 



I fj (p) J d[A] exp str (XA*A - A* AAgfj 



r p 

= J d[g]I p {g) J] sdet-^ /2 (X - gA 3 ) , 



(23) 



where in the last step we performed the integrals over A. Eq. (23) is the 
desired superspace representation of the generating function. Originally, the 
generating function was an integral over ordinary p x n matrices W. The 
representation (23) is an integral over supermatrices g of dimension 4x4 for 
(3 = 1 and 2 x 2 for (3 = 2. This drastically reduces the number of integrals 
to be calculated. 

To complete the description of our result (23), we must discuss the func- 
tion I pig) - This is a super- version of what has come to be called the Ingham- 
Siegel integral [14,15]. (The name is due to Fyodorov [27] who introduced 
such an integral in a related, non-super context.) An important point to ap- 
preciate is that the superdeterminant under the integral sign of (18) depends 
polynomially on the variable <j\. Therefore the Fourier transform (18) does 
not converge in that variable and the supersymmetric Ingham-Siegel integral 
Ip{g) cannot exist as a regular function of g. Nevertheless, we can make sense 
of Ip(g) as a distribution in g, as follows. (In the next subsection, we will 
reproduce the same result (23) by employing convergent integrals only.) 

The Ingham-Siegel integral is invariant under conjugation of a by ele- 
ments of the supergroups UOSp(2|2) or U(l|l) for (3 = 1,2, respectively. 
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This is not automatic, but does hold true once the integral is properly regu- 
larized by an invariant cutoff function. Thus, choosing an invariant cutoff of 
Gaussian form, we define the supersymmetric Ingham-Siegel integral by 



Ip(g):= lim + J d[a] sdet ™^ /2 (l 4 / /3 + iff) exp 



y strap j c (24) 



Here we see the reason why the block 01I2 is multiplied by i = \f— 1: this fac- 
tor cancels the minus sign from the supertrace, thereby making the Gaussian 
cutoff function decrease with increasing real integration variable <j\ . 

By construction, the invariantly regularized Fourier transform (24) is in- 
variant and thus the distribution Ip(g) depends only on the eigenvalues of 
q = URU^ 1 = C/diag(i?o, Ri^/pW^ 1 ■ We note that the diagonal matrix 
Ro has the dimension 2/(3 while R\ is a scalar. An explicit formula for the 
Ingham-Siegel distribution Ip(g) was derived in Refs. [44,15]. The result is 



Ip(g) = K^O(Ro)det^ +1 ^/ 2 - 1 (R Q ) 



:c:iV l-l atrR \(iJL) SiR,), (25) 



with the constants 



K l = -^— K 2 = T ^- (26) 
(n — 2)1 (n — 1)! 

and the Heaviside distribution 

I 1 if R is a positive definite matrix, 
0(i?o) = \o else. (27) 

This completes our description of the result (23). 



3.4 Approach using superbosonization 

We now rederive the result (23) (or, more precisely, an equivalent formula) 
by a different approach, avoiding the use of super-distributions. This will be 
mathematically clean in every respect. The price to pay is that we rely on a 
"black box", namely the superbosonization formula proved in [47]. 

We start over from the very beginning, Eq. (41). Motivated by the fact 
that our probability measure Pp(W, A)d[W] for W G C pxn (resp. W G R pxn ) 
is right-invariant but not left-invariant, we first pass from determinants on 
the left space C p (M. p ) to determinants on the right space C™ (K"), and only 
afterwards introduce the standard Gaussian integral representation: 

(x \ P ~ n det {xil p - WW^) _ det ( Xl l n - W^W) 
\~x~x) det (x l p - WW**) ~ det {x l n - W^W) 
= sdet"' 372 (X ® 1„ - i 4//3 ® W^W) 

= j d[A] exp str (XA^A - AWWAU . (28) 
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X is the same diagonal matrix as before, and we still take Imxo > in order 
for the Gaussian integral to converge. The rectangular supermatrix A is the 
same as before except for the change in dimension p — > n. 

Next we take the expectation with respect to the probability measure 
Pp(W, A)d[W] of the VF-dependent factor under the integral sign: 

J d[W]P {W,A) exp(-jtiWAA^W^ 
= det"^ /2 (l„ ® l p + i/L4 f ® A) 

= sdet~ l3/2 (t 4/l3 (g>l p + iA^A(g> A) . (29) 

In the last step we invoked the duality of Sec. 3.2. For the generating function 
(9) we then get the formula 

Z p {x ,x x ) = (xo/x^-p J d[A]exp ^sbr XA^a\ 

v 

x Y[ sdet^ 72 (U/p + iA^AAj) . (30) 

i=i 

Now a beautiful feature of the integrand in (30) is that it depends only 
on the product A^A which is invariant under left translations An- UA by 
U E O(n) (resp. U(n)) for j3 = 1 (f3 — 2). This means that we are exactly in 
the situation where the superbosonization formula of [16, 17] applies. 

In a nutshell, superbosonization lets us switch from our A-integral to the 
same integral over a supermatrix Q replacing A^ A. The result is 

Zp(x Q ,xx) = {xo/x^-p J DQFp{Q), (31) 

F (Q) = sdet n/3/2 (Q)exp (^strXCjj JJsdet"^ 2 (U/p + iQA,) . 

Q is formally identical to the supermatrix g of Sec. 3.3. The only difference 
is in the domain of integration: the present Q-integral is over a Riemannian 
symmetric superspace (of Cartan type A) with invariant (or Berezin-Haar) 
measure DQ. This is to say that Q runs over the positive matrices in the 
left upper (or boson-boson) block and the unitary matrices in the right lower 
(or fermion-fermion) block. For (3 = 2 the measure DQ is the flat measure 
d[Q] = d[a] of (22) (up to a normalization factor), for j3 = 1 it is the flat one 
times an extra factor of sdet _1 ^ 2 (Q). 

Our two expressions (23) and (31) become identical upon making the 
substitution Q — > \Q/X. In fact, the Heaviside distribution in Rq of (25) 
restricts the integration over all Hcrmitian matrices Qq to the positive ones, 
and the derivatives of the (5-distribution in Ri of (25) have exactly the same 
effect as doing the integral over the unitary variable in the fermion-fermion 
block of Q by the residue theorem. Thus the g-integral in (23) is effectively 
over the Riemannian symmetric superspace parameterized by Q; see [18] for 
more discussion of how the two approaches are related. 
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4 Explicit expressions for the one-point function 

We consider in Sec. 4.1 the complex case (j3 = 2) and rederive the result found 
in Ref. [8]. In Sec. 4.2 we study the real case (/3 = 1) and derive an expression 
as a twofold integral. In both cases we have a choice between the methods of 
superbosonization and generalized Hubbard-Stratonovich transformation. 



4.1 Complex case 

We show how to reproduce in a few more steps the known result [8] for the 
one-point function for = 2. The supermatrix Q in this case is of size 2x2. 
We begin by introducing its eigen-representation: 



Q = URU~ 



U 



— aa —a 
a* \/l — &*a 



R 



r 
s 



(32) 



with eigenvalues r,s and two Grassmann variables a, a*. This parametriza- 
tion is singular (more precisely, degenerate in the Grassmann variables) at 
?' = s. To suppress any effects of the singularity, we are going to utilize the 
scale invariance of DQ to change the integration radius \s\ — 1 to \s\ — q with 
q — > oo. By the change-of-variables formula for superintegrals, our invariant 
integral (31) for the choice of parametrization (32) then takes the form 



J DQF 2 (Q) = (2m)- 1 lim J dr I ds (r - s)- 2 d a d a *F 2 



(33) 



\s\=q 



provided that ^(Q) vanishes on the locus of the coordinate singularity r = 
s = q oo. (Otherwise, so-called Efetov-Wegner boundary terms appear.) 
To arrange for this vanishing property to hold, notice that for large values 
of r = s = q the function -F^glh) behaves as e lc i( x o- x i) _ jf l m (a;o — x\) > 
then this exponential factor makes ^2(^12) vanish on the singular locus for 
q — > +00. We therefore assume the inequality Imio > Imxi. 

Since the two superdeterminants in i^CQ) are functions of the eigenvalues 
r and s only, all of the dependence on a, a* resides in the factor e lstr x ® and 
the process of integrating (actually, differentiating) w.r.t. the Grassmann 
variables is simply done by 



d a d a * e 



istr XQ gi(xor— Xis] 



i(x - xi)(r - s). 



By using (33) and (34) in (31) we obtain 



Z 2 (x ,xi) = (x /xi) n p lim 

q— >oo 



where qa{x\ s) is the function 



1 

Jdrj 



1*1=9 



ds_ r n (x - xi)g A (xi;s) 
2tt s n (r - s) g A (x ;r) 



(34) 



(35) 



g A {x;s) =e- ixs X\{l + isA 3 ). 



(36) 
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Our next step is to perform the inner integral over the circle variable s. 
This is done by invoking Cauchy's integral theorem to show that for any 
complex-analytic function f(s) one has 

J_ I ds = f-r- n f l - n] (r), 1>r, 

27ri J s n {r-s) \+r- n /[ <n l(r), q < r, K > 

where /'-"'(r) stands for the Taylor series of f(r) in r with all terms up 
to order n — 1 deleted, and /^"'(r) = f(r) — /'-"'(r). Application of this 
formula to (35) yields the result 

Z 2 (x ,x 1 ) = -i(x -x 1 )(xo/x 1 ) n -^> ldr 9X } Xl]r } . (38) 

o 

Recall that Imx > Imxi is required in order for this to hold. The correct 
expression for the opposite case of Imxo < Imxi is obtained by replacing in 
(38) every occurrence of i by — i or, equivalently, by sending r — » — r. 

We now turn to the calculation of the one-point function S 2 (x). We recall 
the expression (10) involving Z 2 (x ± ie, X\) and apply our result (38) to the 
case of xq — x ± ie with real parameters x, x\ and e > 0. Our integral 
representations of Z 2 (x ± ie,xi) then combine to a single integral: 



lim \Z 2 (x — ie, x\) — Z 2 (x + ie, ii)) 

e->0+ v y 



dr% -e\r\ 9 A K l2) 4 ( 39 ) 

— CO 

If we re-express the numerator of the integrand as g][ = 9A~ 9a then the 
term g\ contributes (x — x\)8{x — x\) — 0. Hence we replace this numerator 
by —g^ n \xi \ r), which is a polynomial in r of degree n — 1. 

Now for x < we may close the integration contour around the lower half 
of the complex r-plane. From (36) one sees that the integrand is holomorphic 
in r for Imr < 0. We therefore get S 2 (x) = for x < 0. On the other hand, 
for x > the contour has to be closed around the upper half-plane. Again, 
recall Eq. (10). Because the r-integral is now manifestly finite for all values 
of X\ 1 we may safely take the x\ -derivative at X\ = x by simply removing 
the prefactor x\ — x and replacing x\ by x under the integral sign. (Please 
be warned that in (38) the r-integral has to diverge at X\ — x to arrange 
for Z 2 (xq, xq) — 1 in spite of the factor xq — x±.) The result for x > is 



, (40) 



2ttp7 gA{x;r) p 9A(x;r) 



] — >i/Aj 



where the integration contour in the first expression encloses the poles of 
l/gA(x;r) at r — i/Aj (j = l,...,p). To get the second expression, we 

switched from g^ 7 ^ back to — g^ and applied the residue theorem. 
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The residues are best computed by re-inserting the integral representation 
(37) for g^ n \ which gives 

i/4) = (i/A)"- 1 / ~ *- n e- iX3 II (! + ^)- ( 41 ) 



We now use the expansion 



11(1 + isA t ) = ^(ia)*" 1 ^-!^) (42) 

where are the elementary symmetric functions 

E k {A):= A n A 32 ---A 3k . (43) 

i<3l<32<---<ik<P 

The notation A J means that the eigenvalue Aj is to be dropped from the 
diagonal matrix A. Now, inserting the expansion (42) into (41) we carry out 
the s-integral. Our final result for the one-point function then reads 



rt / \ &(x)^ e-^(l/yl,)™ x n ~ k _ . , ? , 

SaW - v g rwi - S g^*-'-*'- (44) 



To conclude this subsection, we present an alternative expression for the 
result (44) as a ratio of determinants: 



S 2 (x) = -^^det 



B A {x) 
C{x) D A 



/AetD A , (45) 
ants {Da) 

Ba{x) and C(x) are row and column vectors with entries 



where Da is the p x p matrix with matrix elements {DA)k.j = A - fe+1 , and 



T n—k 

B A {x) J = e-'^(l/Ajr, C{x) k = -j^-y (46) 

The denominator det Da — Ylj>j'{ A J 1 ~ A J' 1 ) is essentially the Vandermonde 
determinant associated with the numbers A\, . . . ,A„. 

To verify the expression (45) one expands the {p+l) x {p+l) determinant 
with respect to the first row and first column and then uses a standard 
identity [29] for the elementary symmetric functions: 

detgfj = {-ly-'E^M 3 ) (47] 
det D A Ui^O- ~ A ll A i) ' 

where D * is Da with the k th row and j th column removed. In this way 
one immediately retrieves (44) from (45). 

The expression (45) may be useful for certain applications. Also, it is 
easily seen to be directly equivalent to the expression given in [8]. 
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4.2 Real case 

4-2.1 Approach using the eigenvalues of Q 

The main trick in deriving the explicit result for j3 = 2 was to use the eigen- 
representation (32) for the supermatrix Q. We are now going to carry out an 
analogous derivation for (3 = 1. The outcome will be somewhat different in 
that the integrand is no longer meromorphic but has square-root singularities 
in the radial variables of the boson-boson block. 
We use the parametrization 



= k 



Qo " 




k = 


~k 0" 




y/1 — aa 


st 2 _ 




1 2 




a 



fc = 



COS ( 

sin ( 



— sm< 
cos 6 



V V 



(48) 



Qo = 



r a 
n 



which has coordinate singularities at r a = s and r\, = s. As before, we will 
suppress their effects by using the scale invariance of DQ to place the singular 
locus of the coordinate system on the zero locus of the integrand. 

In the coordinates given by (48) our invariant integral (31) takes the form 



/ 



DQFi(Q) = lim 

q— 5-oo 



dr a dr b 



ds 



[0,?) 



1*1=9 
2tt 



167ri (r„ - s) 2 (r b - s) 2 



(49) 



a r b J 27r 



where F\ is the function 

Fi(Q) = sdet n/2 (Q)exp f-strXQ) J^[sdet" 1/2 (1 4 + iQA 3 



(50) 



The expression (49) holds if Fi(Q) and its derivatives vanish on the singular 
locus r a = s = q — > oo and r b = s = q —} co. Observing that our integrand 
Fi(Q) for r a / b = s = q contains the exponential ^i\2 x °~ x ^) 5 we see that F\ 
has the required property if the parameter range is restricted by |lmxo > 
Imxi. We thus impose this restriction. 

It is clear that the integrand (50) does not depend on the angle <j>, and 
all of the dependence on the Grassmann variables is in the factor e lstrX Q/ 2 _ 
Hence the result of doing the Grassmann integral is given by 

d^dndr,, (1 + + V*V) e istrXQ / 2 = e ^(r a +r b )/2-ix lS 
x (i(afo - xi)(r a + r b - 2s) - (x - xi) 2 {r a - s)(r b - s)) . (51) 
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A slight reorganization of the integrand then leads to the following expression 
for the generating function: 

Z lM = (^J% JJ d r ad r bK -n\ / ^ {^y 1 (52) 

[0,9) 2 \s\=q 

Kgo - xi){r a + r b - 2s) _ (x - Xi) 2 \ g A (xi,s) 

(r a - s) 2 {n - s) 2 (r„ - s)(r& - s) j ^/ 5/1 (x ; r a ) 5/1 (x ; r b ) ' 

where g A was defined in (36). 

As before, it is possible to carry out the complex contour integral over 
the variable s. By making partial fraction decompositions to express the 
integrand in (52) by isolated poles (r a — s)~ l and (r b ~s)~ e of degrees £=1,2, 
and then applying the identity (37), we obtain 

7i \ 1 f x A n P ff ^ a Va-r h \^¥~F h n ^ 
Zi(xq,xi) = - — // dr a dr b 



8 \xi) jfjf \/gA{xa;r a )g A (x ;r b ) 

2 

^2(x - xi) v Fy{xi;r a ,r b ), (53) 



+ 

2 

X 

v=l 



where we have taken the upper limits of the r-integrals to infinity, and 

F A V) {x-r a ,r b )= G { A\x;r a )~G^(x;r b ) 

r a - n 

G^\ X ;r)=i^Gf{x;r), G™(x; r) = - 9K - J? T > . (55) 

The result (53) is the analog for ft = 1 of the earlier formula (38) for (3 = 2. 

By recalling the formula (10) we can again deduce an expression for the 
one-point function, S\(x). (This requires a process of analytic continuation 
to remove the restriction Imxo > 2Ima;i.) The result involves only the v = 1 
term from (53) and holds for all n > p+3. Unfortunately, we do not know how 
to do the final two integrals over r a and r b in closed form, and an attempt to 
compute these integrals numerically met with the following complications. 

The first complication is that we do not know how to pass from the ex- 
pression (53-55) involving g > to an analog of Eq. (40) involving g < . Working 
directly with g > is numerically expensive because many terms in the Taylor 
series have to be summed when r a or r b are large. On the other hand, if 
we use the identity g > = g — g < then we incur cancelations, with ensuing 
rounding errors, due to taking the difference of two large numbers. 

The second difficulty is that the exponential part of the integrand os- 
cillates. For the most part, these oscillations can be cured by rotating the 
integration contours for r a , r b to the positive imaginary axis. Ultimately, 
however, these variables still have to go to infinity in the direction of the 
real axis to retain convergence of the integral. The contribution from this 
ultimate part of the integration contour converges rather slowly. 
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These problems notwithstanding, we are still able to verify our formula 
for Si(x). However, by numerical integration using Wolfram-Mathematical 
we are not able to produce stable results in a wide range of parameters. We 
therefore refrain from even writing down that formula, and abandon now the 
eigen-representation (48). Instead, we pursue another approach, exploiting 
the original coordinate system for Q described in Sees. 3.3 and 3.4. 



4-2.2 Direct approach 



Below, we give two further expressions for the generating function. Their 
logical order of presentation depends on which of our two formalisms is used. 

Adopting the generalized Hubbard-Stratonovich transformation, we start 
from the integral representation (23) and evaluate the delta-distributions 
featured in (24-25). We then carry out the Grassmann integrals (or rather, 
derivatives) according to the flat measure specified in (22) to obtain 

1 ff \r a - r b \ ^/r^r b n ~ 3 e~^ a+r ^ /2 

Zi(x a ,xx) = - dr a dr b " x 



llLi V ( x o ~ A l r a )(x Q - Ain) 

l-\) k x p r k ( ^ A 2 ,AfE k JA^)r a r b 

X E V A. [n(n-l)E k (A)+y . 3 ' , 2 . 1 } —. . (56) 
(n-k)\ y j£ (x - Ajr a )(x - Am) 

p / 



3 = 1 

where E k , defined in (43), is understood to vanish when k < 0. The notation 

A jl means that both Aj and Ai are removed from the set {^tz}i=i, 

Alternatively, we can take the formula (31) from superbosonization and 
make the substitution Q i— > iQ/X. Then, by carrying out the Grassmann 
integrals in the original parametrization specified in Sec. 3.4 we obtain 

7 , x (xi/x ) p ff \r a - r b \ y/u¥b n "-' 

Zi(x ,xi) = — -— — dr a dr; 



x </) ds 




167n J J vW^o! ir a /x )g A (xo; ir b /x ) 

g A (xi;is/xi) ( ^ (A J r a )(A j s)(Air b )(Ais) 



A 3 r a )(xi - Ajs)(x - Am)(x 



\ - AjS ( Ajr a A 3 r b 




| s |=i 

-(n-l)V ( -:-- ^— )+,,{„ 

Ajs \x - Ajr a x - Ajr b / 

If we go on to carry out the s-integral, we again arrive at (56) above. Note that 
(56) and (57) make immediate sense as convergent integrals for all X\ € C 
and Xo €E C \ R+. A non-trivial check is Z 1 (x , Xo) = 1. 

From the expression (57) it is easy to verify the large-n behavior pre- 
dicted by the central limit theorem (CLT). Indeed, for large n our integrand 
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develops a saddle point at r a = r& = s —t n due to the presence of the fac- 
tor (y / r a Xb/s)™e~( ra+r ' i> )/ 2+;5 . By steepest descent evaluation of the integral 
around this saddle point one finds 

P x — A 

lim Zi{nxQ,nx\) — I I — , (58) 

n-^roo ±A - Xr\ — Aj 

1=1 ' 

which is the result expected from CLT for the random matrix (WW T )ij — 
Sfe=i WikWjk with independent identically distributed matrix elements Wjk- 



5 Numerical computation of Si(x) 

On recalling the basic identity (10), one immediately produces an expression 
for the one-point function Si (a;); see below. The result looks more involved 
than the one obtained from the eigen-representation of Q, and the reciprocal 
square roots visible in Eqs. (56, 57) remain serious obstacles for further ana- 
lytical calculations. From a numerical perspective, however, the new formula 
has two clear advantages over the previous one: its integrand is free of oscil- 
latory exponentials, and high powers of the integration variables r a and rf, 
occur only through the factor e~( ra+rb ^ 2 , which is numerically inexpensive to 
compute. Thus numerical evaluation of the integral (61) for Si(x) is straight- 
forward, provided that the singularities at r a ,rt, = x/Aj (j = l,...,p) are 
treated with care. We now discuss how to organize this computation. 

To prepare our final formula for Si (x) , we observe that the discontinuity 
Z(x — ie, Xi) — Z(x + ie, Xi) for e — > arises solely from the product of factors 

\/^{x~ A ira ){x- A in ), (59) 

since all other terms in (56) are single-valued across the real x-axis. By the 
derivation of the result (56), the reciprocal square root (59) has the real value 
(x — Air)^ 1 for r a = rt, = r. Its values away from the diagonal r a = are 
defined by the process of analytic continuation. 
For definiteness, let the /1-values be ordered by 

< Ai < A 2 < . . . < A p . 

Fixing any x > 0, we partition (a dense open subset of) the domain K + for 
r a and into p + 1 intervals Uq(x) := (0,x/A p ), U p (x) := (x/ A%,oo), and 

U p ^{x):^{x/A l+U x/Ai) (l = l,...,p-l). (60) 

The product over i = 1, . . . ,p of the square roots in (59) is then discontinuous 
across the s-axis only for (r a , rj,) S Ui(x) x Uy [x) with l + l' an odd number. 
Thus we arrive at an integral representation for Si(x) of the form 

dr b / ' (61) 

Y[U^-A t r a ){x- A in )\ 



Si(x) = j dr a ji 

Ui(x) U t , 
fx,A( r a,r b ) + 



- Ajr a ^T 1 ., (x - A J r a ){x - A j ,r b ) i 



j'=i i^j 
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The outer sum is over the pairs (/, I') £ {0, 1, . . . ,p} 2 with odd sum I + V , 
and the fx\{ r at r b) are certain real-valued analytic functions, which do not 
need to be written down as they are easy to read off from (56) using (10). 

What we do need is some further discussion of the precise meaning of the 
integral (61). There is no problem with the first summand (containing f x ,A), 
as the square root singularities from (59) are integrable. However, in the 
second and third summand of (61) the order of the singularity is enhanced 
from (x — Ajr)^ 1 / 2 to (x — Ajr)~ 3 / 2 for one (j) or two (j ^ j') factors of 
the product over i = 1, .. . ,p. These terms have to be properly understood 
as principal-value integrals by their definition as a limit e — > 0. To render 
them suitable for numerical evaluation, we use a partial integration identity 
of the sort 



— lim 

2 £->o 



g(r) dr 
(r - ie) 3 / 2 



' 7-1/2 



+ 



r— — 5- 



d r g(r) dr 

r !/2 : 



(62) 



which holds for any smooth function g(r). We apply this identity to the 
present situation with r = r a — Aj/x and/or r = rj, — Ay jx (and suitable 
values for 5±), and thereby regularize the integral around the 3/2 singulari- 
ties. In order for this rcgularization to be well-defined, we need the diagonal 
Ta = fb to lie (as it does) outside the domain of integration, so that our 
integrand g(r) containing the factor \r a — rf,| is in fact smooth. 
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Fig. 1 One-point function Si (x) for the parameter values p = 10 and n = 50. The 
set of {yli} ; - = i,...,io is {1,0.49,0.4225,0.36,0.25,0.09,0.0729,0.0529,0.04,0.0225}. 
The solid line is the result obtained by numerical integration of the analytical 
formula (61). The bin size of the histogram is 0.7 for the large figure and 0.1 for 
the inset. The inset magnifies Si(x) near the origin. 
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Fig. 2 Si(x) for the same values of p and {A,} but a larger parameter n = 200. 
The solid line is again the result obtained from Eq. (61). The bin size is 2 for the 
large figure and 0.3 for the inset. Note that because of the bigger value of n the 
peaks are more pronounced (than those in Fig. 1) around the CLT- values nAj. 

By using the commercial software Mathematica® [28] we now numer- 
ically evaluate the integrals in Eq. (61). In Figs. 1 and 2 we show examples 
for (p, n) = (10, 50) and (10, 200), respectively. Our result (solid line) is com- 
pared with a Monte-Carlo simulation (histogram) using an ensemble of 10 5 
random matrices. We see that the agreement between the Monte-Carlo sim- 
ulation and our result from (61) is perfect. 

We observe that the last peak, centered around x = 50 in Fig. 1 and 
x = 200 in Fig. 2, lies near nA p = n, in agreement with what is expected by 
the central limit theorem (CLT). The other peaks are slightly shifted from 
their asymptotic CLT-positions nAj. With increasing n these shifts become 
smaller and the peaks become more pronounced, as it should be. 



6 Summary 

We have developed a supersymmetry approach to derive exact expressions 
for the one-point function of real and complex Wishart correlation matri- 
ces. Supersymmetry got us around a difficult group integral which arises in 
the traditional approach for the real case. The crucial advantage of the su- 
persymmetry approach is the drastic reduction of the number of integrals. 
For both cases, complex and real, we showed how to express the one-point 
function as an eigenvalue integral. In the complex case the supermatrix has 
dimension 2x2 and thus two eigenvalues. We carried out both integrals and 
demonstrated the equivalence of our result to that of Ref. [8] . 

In the more demanding real case (/3 = 1), the supermatrix has dimension 
4x4 and 3 distinct eigenvalues. While the Efetov-Wegner boundary terms 
due to diagonalization for /3 = 1 have never been given in explicit form, 
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we showed how to suppress them by a variable substitution that pushes 
the coordinate singularities outside the domain of integration. One of the 
three eigenvalue integrals is easily done by residue calculus. We do not know 
how to calculate the remaining twofold integral by analytical means, and 
an attempt to compute it numerically met with some complications. We 
therefore abandoned the coordinate system given by the eigenvalues and 
turned to a direct approach using standard coordinates. Thus we produced 
a second formula for Si(x), still as a twofold integral, which we were able 
to compute numerically in a stable and efficient way. We also illustrated our 
result by comparison with a Monte-Carlo simulation. 

Previous approaches to the real case had resulted in slowly converging 
series of Jack or zonal polynomials. We believe that our formula represents 
a considerable improvement over these previous results. From a conceptual 
viewpoint, one might say that our result re-sums a multiple infinite series of 
Jack polynomials (and integrals thereof) in a non-trivial real case. 

Thus we hope that we have demonstrated that the supersymmetry method 
is a powerful tool to tackle problems in multivariate statistics. 
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